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ABSTRACT 


Two main problems are analyzed and discussed in this thesis. First, a 
detailed derivation of the exact solution for the inhomogeneous Helmholtz 
equation for a free-space propagation problem when the square of the index of 
refraction is a linear function of depth and the input is an omnidirectional point 
source is performed and discussed. The exact solution is in terms of Airy and 
Bairy functions. Second, the accuracy of the Recursive Ray Acoustics (RRA) 
Algorithm was tested by generating test cases using a sound-speed profile with the 
square of the index of refraction a linear function of depth. The acoustic pressure 
(amplitude and phase) calculations obtained from the RRA Algorithm were then 
compared to the exact "Airy function" solution. Computer simulation results 


indicate that both solutions are in good agreement. 
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I. INTRODUCTION 

This thesis has two primary goals. The first goal is to derive and document the 
solution to the three-dimensional inhomogeneous Helmholtz equation for a free-space 
propagation problem when the square of the index of refraction is a linear function of 
depth. Given this type of ocean medium, the inhomogeneous Helmholtz equation has an 
exact solution in terms of Airy functions. This exact "Airy function” solution will be 
incorporated into the computer program Linear Space-Variant Ocean (LSVOCN) as an 
additional ocean-medium transfer function. The second goal of this thesis is to further 
test the Recursive Ray Acoustics (RRA) Algorithm by comparing the magnitude and 
phase of the acoustic pressure calculated along a ray path by the RRA Algorithm with the 
magnitude and phase calculated by LSVOCN. 

The first step in this process requires a careful derivation of the solution to the 
three-dimensional inhomogeneous Helmholtz equation for a free-space propagation 
problem when the square of the index of refraction is a linear function of depth. This 
derivation will be discussed in detail in Chapter II which is divided into three main 
sections. Section A defines the three-dimensional inhomogeneous Helmholtz equation to 
be solved. (By using a zeroth-order Hankel transform, the solution process is simplified 
to solving a one-dimensional inhomogeneous ordinary differential equation vice a 
three-dimensional inhomogeneous Helmholtz equation). Section B defines the 
mathematical formula used to model the square of the index of refraction as a linear 
function of depth. Section C documents the three steps that lead to the final solution of 
the inhomogeneous Helmholtz equation. 

First, the description of the sound-speed profile is defined. This description 
allows us to solve for the two constants in the equation for the square of the index of 
refraction as a linear function of depth. Second, the homogeneous depth equation is 
manipulated into the form of Airy's differential equation and then solved in terms of Airy 


and Bairy functions. Third, the solution of the inhomogeneous depth equation is obtained 


by applying the appropriate boundary conditions at the source. 





The solution to the inhomogeneous Helmholtz equation is divided into two cases. 
In case one, the square of the index of refraction has a negative gradient. This case is 
further subdivided into a solution for the receiver depth equal to or above the source 
depth, and a solution for the receiver depth below the source depth. In case two, the 
square of the index of refraction has a positive gradient, and again a solution for the 
receiver depth equal to or above the source depth, and a solution for the receiver depth 
below the source depth are found. 

Chapter III is divided into three sections. Section A discusses the numerical 
techniques used in the process of converting the Airy function solution into working 
FORTRAN computer code used by LSVOCN, and defines the Green's function used to 
test LSVOCN. Before testing of the RRA Algorithm began, a comparison was made of 
the output of LSVOCN, using the above derived Airy function solution, with a test case, 
in an effort to validate the FORTRAN computer code. The simplest test case considered 
compared the magnitude and phase of the time-independent free-space Green's function 
for an isospeed medium with the output from LSVOCN. 

Section B of Chapter III discusses the three different Green's function test cases 
that were used to test LSVOCN. These are: (1) Green's function test case with the source 
sound-speed varied such that the modeled ocean medium approached an isospeed 
medium in order to compare the magnitude and phase of the Green's function with 
LSVOCN, (2) Green's function test case with the horizontal range varied while 
maintaining the receiver at the source depth, thus simulating an isospeed ocean medium 
for comparison of the Green's function magnitude and phase with LSVOCN, and (3) 
Green's function test case with the receiver depth varied above and below the source 
depth by one-half meter allowing for an approximate comparison of the Green's function 
magnitude and phase with LSVOCN. In each test case, the change in sound speed 
between the source and receiver was small enough to approximate an isospeed ocean 
medium. 

Section C of Chapter III focuses on the validation of the predicted values for the 
magnitude and phase of the acoustic pressure along a ray path given by the RRA 





Algorithm. This is done by propagating sound rays at various launch angles using 
negative and positive gradients for the index of refraction profile. The horizontal range 
and depth coordinates of selected points along these ray paths have been recorded and 
will become the receiver horizontal range and depth coordinates used later as inputs to 
LSVOCN. The magnitude and phase values calculated by LSVOCN are then compared 
to the corresponding magnitude and phase values obtained from the RRA Algorithm for 
validation purposes. 

Chapter IV summarizes the numerical problems encountered during the validation 


process, states conclusions as well as recommendations for future research. 





Hf. AIRY FUNCTION SOLUTION OF THE INHOMOGENEOUS 
HELMHOLTZ EQUATION 


A. THE INHOMOGENEOUS HELMHOLTZ WAVE EQUATION 
In this chapter we shall derive a solution for the following inhomogeneous 


Helmholtz equation in cylindrical coordinates (r, @ ,y), shown in Fig. 1: 
2 2 5 r 
Fone.)+ 72.07 (ry) + Z0dr9) + POAT) =F28Y-yo) (21) 


where the velocity potential p,(r,y) is axisymmetric (i.e., independent of the azimuthal 


angle @ ), 


Ky) = (2.2) 


is the depth-dependent wave number, c(y) is the depth-dependent speed of sound, and the 
impulse functions on the right-hand side of Eq.(2.1) represent a unit amplitude, 


omnidirectional point source at r= 0 and y= yp. 


(r,9,y) 





Figure 1. Illustration of problem geometry. 








The solution 9,(r,y) of Eq.(2.1) can be expressed as an inverse zeroth-order 


Hankel transform as follows, see [Ref. 2, Appendix 3C]: 


OAT,Y) = Ho {Pkr,Y)} =] Ok, yolk rkrdk, (2.3) 
where 
Ok-,Y) = HolOdr.y)} = |” or, yolker)rdr (2.4) 


is the forward zeroth-order Hankel transform (also known as the Fourier-Bessel 
transform) of @,(7,y). Therefore, if we can find ®,(k, y), then substituting ®,(k,,y) into 
Eq.(2.3) will yield a solution of Eq.(2.1). 


Since 
io Z04r.y) + 12047») =-BOMk,,») (2.5) 
or? > r or > r ro . 
and 
Ho( 22) a1, (2.6) 


taking the zeroth-order Hankel transform of Eq.(2.1) yields 


oy — 
rm 5 Ar,3) + KO kr, y) = APO) a (2.7) 
where 
KV) =?) - 7. (2.8) 


In Eq.(2.8), k,(y) is the depth-dependent propagation-vector component in the Y 
direction and k,is the constant propagation-vector component in the horizontal range r 
direction (see Fig. 2). Therefore, now we only have to solve the ordinary differential 
equation given by Eq.(2.7) instead of the partial differential equation given by Eq.(2.1). 
Substituting the solution ®,(k, ,y) into Eq.(2.3) yields the solution Q,(r,y) to Eq.(2.1). 





Figure 2. Illustration of the propagation-vector components k, and k, (y) 
and their relationship to the wave number f(y). 


B. PROBLEM CONSTRAINTS 
The solution of Eq.(2.1) that shall be derived is only applicable to free-space 
propagation when the square of the index of refraction is a linear function of depth y. The 


index of refraction is defined as follows: 


ny) = =o (2.9) 


where c, = c(y,) is the speed of sound (m/s) at the source depth y, and c(y) is the speed 
of sound (m/s) at depth y. If we let k, = k(y,), then evaluating Eq.(2.2) at y= y, yields 


ko =f (2.10) 


Therefore, we can write that 


Ky) = konty) =. (2.11) 


As a result, Eq.(2.8) can be rewritten as 








2 (y) = kon?(y) -K2. (2.12) 
To ensure that the square of the index of refraction is a linear function of depth, we 


define the following equation: 
n*(y)=aiy+ao. (2.13) 


C. PROBLEM SOLUTION 

Using the results from the last two sections, specifically Eq. (2.7) and (2.12), we 
are now ready to solve the inhomogeneous Helmholtz equation. The solution can be 
broken into two parts. First, determine the solution of Eq.(2.7) at y= yy, i.e., solve the 
homogeneous equation. Solving the homogeneous equation will result in four unknowns. 
In order to determine these four unknowns, which is the second part of the problem, we 
apply the boundary conditions above and below the source and the boundary conditions 
at the source, i.e., at y = y,. As shown in [Ref. 2, Appendix 3C], the two boundary 
conditions that must be satisfied at the source are continuity of acoustic pressure and 
discontinuity in the y component of the acoustic fluid velocity vector. With this brief 
introduction in mind, let us now proceed with the solution of the inhomogeneous 


Helmholtz equation. 
1. Description of the Sound-Speed Profile 


Consider a sound-speed profile where the square of the index of refraction is a 


linear function of depth, [see Eq.(2.13)]: 
n*(y) =ayy+ao. (2.14) 


With the above equation in mind, we must now solve for the two unknown constants, 


a, and a,, in order to solve Eq(2.7). Substituting Eq.(2.9) into Eq.(2.14) yields 





co =ayt+a (2.15) 
Ghee 





Rewriting Eq.(2.15) yields 








1 a ao 
= Ay 4 2. 2.16 
AQ) ote ie, 
and if we let 
a, = 41 (2.17) 
Co 
and 
al, =“, (2.18) 
Co 
then 
1 / 
AG) =aiyt+da. (2.19) 


Solving for the speed of sound yields 


(2.20) 


1 Co 
cy) = ———. 
{aiytag {1+ 


Equation (2.20) will allow the sound-speed profile to be calculated, once the unknown 


constants have been determined. 


Evaluating Eq.(2.19) at y = 0 and at y = y, , i.e., the source depth, yields 








Thus, using Eqs.(2.21), (2.23), (2.17), and (2.18), a, and a, can be determined. 

For our problem, we will use a source depth of y, = 150 meters. In order to 
define a sound-speed profile, we specify the sound speed (m/s) at y = 0 to be c(0) = 1500 
m/s. For a positive sound-speed gradient we specify the sound speed at the source depth 
to be c, = 1515 m/s. Similarly, for a negative sound-speed gradient, we specify the 
sound speed at the source depth to be c)= 1485 m/s. These sound-speed profiles, along 
with the corresponding index of refraction plots, are shown in Fig. 3 and Fig. 4, 


respectively. 


INDEX OF REFRACTION SQUARED vs. PEP SOUND SPEED PROFILE 
0 


00 0 L i 
3°) 96 1 702 {6004610 46201690 
INDEX OF REFRACTION SQUARED & (rs) 





Fig. 3 Positive sound-speed gradient. 
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Fig. 4 Negative sound-speed gradient. 
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2. Solution to the Homogeneous Depth Equation 


Having carefully defined the sound-speed profile, we are now ready to solve the 
inhomogeneous depth dependent Helmholtz equation as given by Eq. (2.7). As stated 
previously, the first part of the solution process is to solve the homogeneous equation 
where the square of the index of refraction is a linear function of depth. The 


homogeneous equation is 


Fx OAky) + BOVO(Ey) =0. (2.24) 


Substituting Eq. (2.13) into Eq. (2.12) yields 
1) = Blay+ao)—#, | (2.25) 


or 
K2(y) = kKoary + kao — k?. (2.26) 


Substituting Eq. (2.26) into Eq. (2.24) and simplifying 
Sx Ok) + (arkay + a0k} — YO (k,,y) =0. (2.27) 


If we let & and &, be defined as 


Oto = aok? — k? (2.28) 
and : 

O11 = a1 ke, (2.29) 
then 

Fz Pky) + (Oy + Co) Pky) = 0. (2.30) 


Next, we must manipulate Eq. (2.30) into the form of Airy's differential equation. 





The first step in this process is to define the new term 


2 
Ci) =+(a1) 3 (Oy + Qo). 
Applying the chain rule yields 


d aC(y) 
works, Nz nt me 


and since 


upon substituting Eq.(2.33) into Eq. (2.32), we obtain 





4 ok, y) = #(011)3 2-0 (ky). 
dy dv) 


Applying the chain rule again yields 


£ OKk,,y) = 41 (k,,y) = 4 sta )3 


dy? Be dy 





=+(04)3 LoKk,,y), 


EO dy 


and upon substituting Eq. (2.34) into Eq. (2.36), we obtain 





Z 2 
yr elkny) =+(a1)3 O/k,,y). 


d? 
a? (y) 
Combing Eggs. (2.30), (2.31), and (2.37) yields 


ry) + (4)(01)3 CO)OAk,,¥) = 


Ho4)3 


5) de7gy 








d 
O/k,, 
at(yy | 


0, 


(2.31) 


(2.32) 


(2.33) 


(2.34) 


(2.35) 


(2.36) 


(2.37) 


(2.38) 


or 


Ar gy ier £9) =0. (2.39) 


Upon choosing the minus sign, we have Airy's differential equation: 


Fite) Cy) (k-, y) = 0. (2.40) 


The general solution of Airy's differential equation is given by 


Ok,,y) = AAi[C(y)] + BBiI[G()] , for y < yo (2.41) 


®(k,,y) = CAi[C(y)] + DBi[C(y)] , for y > yo (2.42) 


where Ai and Bi are the Airy and Bairy functions, and A, B, C, and Dare unknown 
constants. Equation (2.41) and Eq.(2.42) form the solution to the depth dependent 
homogeneous Helmholtz equation. We have an equation with four unknown constants as 
a result of solving a second-order differential equation. To find these four unknown 
constants, we must apply the boundary conditions above and below the source and the 
boundary conditions at the source, as will be done in the next section. 

3. Solution to the Inhomogeneous Depth Equation 

In this section we will determine the general solution to the inhomogeneous depth 
dependent Helmholtz equation. As stated previously, in order to solve for the four 
unknown constants in Eq.(2.41) and Eq.(2.42), we must satisfy the boundary conditions 
above and below the source and the boundary conditions at the source. However, before 
we proceed with this solution, let us examine in detail the behavior of the functions 
Ai{C(y)] and Bi{C(y)] to gain insight into the solution of the problem. 


By choosing the minus sign in Eq.(2.31), we can write that 


_2 
CO) =-{01) 3(a1y +00) (2.43) 





Because 


Ky) = any + ao, (2.44) 
Eq.(2.43) can be rewritten as 
2 
G0) =— (01) 342(y). (2.45) 


Furthermore, by referring to Eq.(2.12), we can write that 
k5(y) = koln?(y) — sin?(Bo)] (2.46) 


since k, = k,sin(B,), where B, is the launch angle, measured with respect to the positive Y 
axis. Values of B, range from 0 to 180 degrees, [Ref. 2, Chapter 5]. Tables 18 through 
23, located in the Appendix, contain values for n(y), k,/(y), and C(y) as a function of 
frequency, B,, and depth y, for the analysis of the functions Ai{E(v)] and Bi[C(y)]. The 
data in these tables is further displayed in graphical form in Figures 12 through 17, also 
located in the Appendix. 

To analyze the behavior of Ai[C(y)] and Bi[C(y)] , let us draw our attention to Fig. 


5 and Fig. 6, as shown below: 





Figure 5. Case one: Figure 6. Case two: 
n’(y) with negative gradient. n(y) with positive gradient. 





The free-space solution can be broken into two cases. Case one: n’ (y) has a negative 
gradient, (see Fig. 5, Fig. 3, Fig. 15 through Fig. 17, and Tables 21 through 23). Case 
two: n’(y) has a positive gradient, (see Fig. 6, Fig. 4, Fig. 12 through Fig. 14, and Tables 
18 through 20). Also, as can be seen from Fig. 5 and Fig. 6 , each case must be analyzed 
in two regions, one above the source depth y, , and one below the source depth yy. 

First, let us analyze the behavior of the functions Ai[C(y)] and Bi[E(y)] for, 
Region IIb, Case one. Notice from Fig. 5, that the value of n’ (y) is less than one in this 
region. Thus, as can be seen from Eq.(2.46), k(y) will take on negative values. This will 
result in values of C(y) taking on positive values [see Eq. (2.45)]. Therefore, in Region 


IIb, Case one we should use the function Ai[C(y)] to ensure a bounded solution since 


Ai[C(y)] ------> 0 as C{y) ------> too 


and 


Bi[(Q)] > +0 a8 C9) ——> +=. 


Following a similar argument for Region ITa, Case two we should only use the function 
Ai[C(y)]. Thus, by discarding the Bairy function in Region IIb, Case one and Region IIa, 
Case two we have ensured that our solution will remain bounded above and below the 
source as C(y) approaches positive infinity. 

Second, for Region Ila, Case one and Region IIb, Case two the value of n’ (y) is 
greater than one. Thus kV) will take on positive values, which leads to negative values 
for C(y). Therefore, both Ai[G(y)] and Bi[C(y)] should be used in Region Ila, Case one 
and Region IIb, Case two since both the Airy and Bairy functions are well behaved for 
negative arguments. The asymptotic forms of the Airy and Bairy functions are defined as 
follows, [Ref. 4, Chapter 7] 


Amt = Fel? sino? +71 LO) > te, (2.47) 





BiLO] = ILO) cosl FLO +2] Ey) —> 40. (2.48) 


Thus, to determine the general solution of the inhomogeneous depth equation we 
must ensure the solution is bounded in Region IIb, Case one and Region IIa, Case two, 
i.e., D = 0 (Case one), see Eq.(2.42) and B = 0 (Case two), see Eq.(2.41). Next, we must 
ensure the solution is in the form of a propagating wave in Region IIa, Case one and 
Region IIb, Case two, i.e., A = B, and a = B/A (Case one), see Eq.(2.41), and C = B, and 
b = D/C (Case two), see Eq.(2.42). In summary, we have the following forms for the 
free-space propagation solutions for each case based on Eq.(2.41) and Eq.(2.42): 


Case one: ®,,(k,.¥) = B, [A C)] + aBi[Co)]] fory < y, (2.49) 
®,.(k, ¥) = A, [ATCO for y 2 yo (2.50) 
Case two: ®,,(k, ¥) = A, [ATS] fory< yy (2.51) 


Py o(k,v) = B, [ATS] + SBTC] for y 2 y, (2.52) 


where A, and B, are arbitrary constants and a and b are known constants to be discussed 
later. Thus, for each case, we have two equations and two unknowns. 

To solve for the two unknown constants, we apply the boundary conditions of 
continuity of acoustic pressure and a discontinuity in the y component of the acoustic 
fluid velocity vector at the source depth y,: 


1. continuity of acoustic pressure: 
Dp5(k-, Yo) i Drak, yo) =0 (2.53) 


2. discontinuity in the y component of the acoustic fluid velocity vector: 


J br(kr Yo) ~ Jp Pralkr Yo) =a (2.54) 





where d = x is the finite discontinuity, [Ref. 2, Appendix 3C]. 


CASE ONE: 
Substituting Eq.(2.49) and Eq.(2.50) into Eq.(2.53) yields 


AyAi[C(vo)] — By[Ai[C(vo)] + aBilC(vo) ]]= 90, 
and upon solving for A, we obtain 


[AilCQo)] + aBiSQo)]] 
AilE(vo)] 


Ay =By 
Substituting Eq.(2.49) and Eq.(2.50) into Eq(2.54) yields 
(0) 5AyAi? [S(v0)] - (-)(011)3B, LAH [(0)] + aBi' [Go =4, 


and upon substituting Eq.(2.56) into Eq.(2.57) and rearranging terms, we obtain 


aBy ‘ -/ _ of = d 
FAG Gon BEOMAMEOON ~ AiCOOIBF LOW =A. 


Since, see [Ref. 1, Chapter 4] 


BilC(yo) AH [G(r0)] ~ AilSro)]BF16V0)] = - 


(2.55) 


(2.56) 


(2.57) 


(2.58) 


(2.59) 


which is known as the Wronskian of the Airy and Bairy functions, substituting Eq.(2.59) 


into Eq.(2.58) and solving for B, yields 





B, = 4 _aiteoo)). 
a(a1)3 





(2.60) 





Substituting Eq.(2.60) into Eq.(2.56): 

nd : ; 
TIAMC(vo)} + aBi[C(vo)]]. (2.61) 

a(o,)3 





ASE TWO: 
Substituting Eq.(2.51) and Eq.(2.52) into Eq.(2.53) yields 
(2.62) 


B,{Ai[C(vo)] + bBi[E(vo)]] — AyAi[C(v0)] = 0, 


and upon solving for A, , we obtain 
[Ail(vo)] + SBilC(yo) I] 
A,=B ; 2.63 
oY AICO) — 
Substituting Eq.(2.51) and Eq.(2.52) into Eq.(2.54) yields 
(2.64) 


~(011)3ByfANIC(v0)] + BB/[C(y0) 1] — (ers) 34,A"{CO0)) =, 





and upon substituting Eq.(2.63) into Eq.(2.64) and rearranging terms, we obtain 
(2.65) 





bBy . / _ R J Bien 
Fira AUB (Loo)) — BilSro)I4F [CoN a 
Since 
Ai{C(vo)IBUC(v0)] - BC) 1Ai Co) = 4, (2.66) 
substituting Eq.(2.66) into Eq.(2.65) and solving for B, yields 
By=-—™4 TAC 0)). (2.67) 
b(a1)3 





Substituting Eq.(2.67) into Eq.(2.63) yields 





Ay=-—™ZJA,[C00)] + SBC 0)I. (2.68) 
b(a1)3 


Values for the constants a and b must be determined before we have a complete 
solution. To determine the value of a, recall Eq.(2.49) and Eq.(2.60). Substituting 
Eq.(2.60) into Eq.(2.49) and rearranging terms yields 


Onakr,y) = AWOL eon] + BACON], vs 9 (2.69) 


(a1)3 


Note that ( see Tables 21 through 23, located in the Appendix) as y > — © , C(y) 4 — 00. 
Recall the asymptotic expansions of Ai{C(y)] and Bi[C(y)]: 


AICO] = Felon sinkg-ton}? +4], (9) —> +20 (2.70) 
and 
BilQ)) = FELON? costdConl? +3, LQ) > 4. (2.71) 


Therefore, substituting Eq.(2.70) and Eq.(2.71) into Eq.(2.69) yields 


Dpelkr,y)= EON ony] EsinGiLonl? +H + cosl-tonl? + »|.2.72) 


J (01)3 


Since this derivation corresponds to a time-harmonic solution with time dependence 
given by exp(+/2nfi), in order to get ®,,(k, ,y) into the form of a wave traveling in the 


positive -C(y) direction, we must choose a = +/ since 


exp(-jI-()]) = cos([-C()}) —/ sin([-0)]). (2.73) 








To determine the value of b, we follow a similar logic. Substituting Eq.(2.67) into 


Eq.(2.52) and rearranging terms yields 


Das(kr,y) = EOI 1 aiteyy] + BiCON], yey. 
3 


(Q;) 


(2.74) 


Note that (see Tables 18 through 20, located in the Appendix) as y+, C(y)  — 0, 
Using the same logic as before, since -C(y) --> +o, we must pick b = +j to ensure the 
correct form of Eq.(2.74). 


and 


and 


In summary, we have the following formulas: 


CASE ONE: 
D pal, y) = OO) sire(y) 1+ jBACONI} .y <¥0 
(47)(011)3 





@polk,,y) = nd{Ai[C(vo)] + JBAG vo} AiICON > yo 
(+7)(01)3 
CASE TWO: 
Opal y) =~ EAAMLOW *SBMCO grey y ye 
(47)(01)3 
Dpn(k,y) = —- AO aiteyy BACorl) .y>ye 
(+7)(Q;)3 
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(2.75) 


(2.76) 


(2.77) 


(2.78) 





Recall that d= sills, 
20 


Now that a solution for ®,(k, ,y) has been derived, all that is left to determine the 
velocity potential is to insert ®,(k, y) into Eq.(2.3) and integrate. It is interesting to 
note, as described in, [Ref. 2, Chapter 9], that: 


Ok, ¥) =H m(FSr.y 5 Yo) (2.79) 


where H_(ff.y ; y, ) is the time-invariant, space-variant ocean medium transfer function. 
Thus, as described in the Introduction of this thesis, Eq.(2.75) through Eq.(2.78) will be 
added to the already existing library of ocean medium transfer functions in the program 
LSVOCN. With this derivation and new transfer function, we can compare the results 
from LSVOCN to the RRA algorithm to determine the accuracy of the RRA algorithm. 








Ill. COMPARISON OF THE "AIRY FUNCTION" SOLUTION 
(MAGNITUDE AND PHASE) 
WITH THE RRA ALGORITHM'S SOLUTION 


A. THE TIME-HARMONIC FREE-SPACE GREEN'S FUNCTION AND 
NUMERICAL ANALYSIS TECHNIQUES 
1. The Time-harmonic Free-Space Green's Function for a Homogeneous 
Ocean Medium 

Once the solution to the inhomogeneous Helmholtz equation for a free-space 
propagation problem with the square of the index of refraction a linear function of depth 
had been derived, the solution was converted to FORTRAN computer code and added 
to the computer program LSVOCN as a new ocean medium transfer function. LSVOCN 
will accept as primary inputs: 

1. Sound-speed at the ocean surface and source, in meters/second. 

2. Receiver depth and horizontal range relative to the source, in meters. 

3. Frequency the source will operate, in Hertz. 

The output of LSVOCN is either the velocity potential (magnitude and phase) or 
the acoustic pressure (magnitude and phase) at a specified receiver location. To test the 
RRA Algorithm, the depth and horizontal range coordinates of selected points along a 
ray path, propagated by the RRA Algorithm, will be input into LSVOCN as receiver 
coordinates. The calculated acoustic pressures (magnitude and phase) from LSVOCN 
corresponding to these selected points, will be compared with the calculated acoustic 
pressures (magnitude and phase) from the RRA Algorithm. 

Before beginning the testing process, we needed to ensure that both the 
theoretical solution and FORTRAN computer code were correct. In order to perform this 
validation, test cases were developed for which theoretical values of the velocity 
potential could be calculated and compared with the predicted values found using 
LSVOCN. The time-harmonic free-space Green's function which corresponds to a 
spherical sound source operating in the monopole mode of vibration in a homogeneous 


ocean medium, was selected as the basis for these test cases. 








As shown in [Ref. 2, Chapter 4], the velocity potential of the acoustic field 
produced by a spherical sound source in the monopole mode of vibration operating in an 


unbounded, homogeneous medium is given by 


exp(+jka) 


> 
ae R>a, (3.1) 


OAR ,6 , W) = Sog/rlo) 


where S, is the source strength (volume flow rate) in m’/s, 


_ _ &xp(-/kR) 
gril) =- a (3.2) 


is the time-harmonic free-space Green's function (exp(jax) factor suppressed) and R is 
the spherical range between source and receiver. In the limit as a > 0 (modeling a point 
source) Eq.(3.1) becomes 


—jkR 
GAR 0 ,W) = SogArl0)=—S PLY | p> 0. (3.3) 


Using Eq.(3.3), theoretical values for the magnitude and phase of an acoustic 
field at a receiver in an unbounded, homogeneous ocean medium located R meters from a 
point source can easily be calculated. Although these values are valid only for the 
homogeneous case, as the modeled ocean medium that LSVOCN uses becomes more 
homogeneous, then the calculated velocity potential (magnitude and phase) from 
LSVOCN should approach the predicted values calculated by Eq.(3.3). The theoretical 
magnitude of the velocity potential calculated by Eq.(3.3) is given by 


_ So. 3.4 


The theoretical phase of the velocity potential calculated by Eq.(3.3) is given by 
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LOr= n( 2 ~ e (3.5) 
where fis the frequency at which the source operates, R is the spherical range between 
the source and receiver, and c is the speed of sound. In order to test the Airy function 
solution and the corresponding FORTRAN computer code in LSVOCN, Eq.(3.4) and 


Eq.(3.5) were used in three different Green's function test cases, defined as follows: 


1. With the receiver and source depth equal, and with the receiver one meter from 
the source, vary the speed of sound at the source so that it approaches the speed 
of sound at the surface. 


2. With the receiver and source depth equal, and the source sound-speed constant, 
vary the horizontal range between the source and receiver. 
3. With the source sound-speed constant, vary the receiver depth. 


Each of these test cases will be covered in detail in Section B of this chapter. 
2. Numerical Analysis Techniques 


Recall from Eq.(2.3) that 
Ofr,y) =H {Odk,,y)} =], Ok, yolk rkrdky. (3.6) 


In order to calculate the exact "Airy function" solution at a given receiver depth and 
horizontal range, the integral in Eq.(3.6) must be computed. As we discovered, this 
integral is difficult to evaluate numerically. The integrand in Eq.(3.6) is oscillatory due 
to the Bessel function, when negative arguments of the Airy and Bairy functions are 
involved. Numerous overflow and underflow errors were generated when Eq.(3.6) was 
calculated numerically with the computer. To eliminate overflow errors, the 
exponentially scaled versions of the Airy and Bairy functions were used. The following 
example illustrates the method of using the exponentially scaled Airy and Bairy functions 
to eliminate overflow errors. | 

In order to solve Eq.(2.75) or Eq.(2.78) for the case of propagating waves, it is 
necessary to multiply the Airy and Bairy functions together. Although the argument of 


these functions, C(y), is more likely to take on negative values in regions Ila, Case one 








and region IIb, Case two, once the propagating waves turn evanescent, the value for C(y) 
is forced to take on positive values which cause the Bairy function to blow up, leading to 
overflow errors. By experimentation it was discovered that the computer could handle 
values of C(y) up to positive forty, without generating overflow errors. However, to 
include most of the evanescent waves, C(y) must be allowed to take on much larger 
positive values, which requires the use of the exponentially scaled Airy and Bairy 


functions: 


1.DAI(x), DBI(x) -- double-precision IMSL Airy and Bairy functions. 


2.DAIE(x), DBIE(x) -- double-precision IMSL exponentially scaled and Bairy 


functions. 
where 
DAIE(x) = DAI(x) x exp(2x2) (3.7) 
and 
DBIE(x) = DBI(x) x exp(— 2,3), (3.8) 


Therefore when multiplying Ai[C(y)] and Bi[C(y)], first form the difference term: 
2 2 
DIFF = al t0)3 —E(y)3 i (3.9) 


Then perform the required multiplication as follows: 


A,[S(0)] x Bi[GQ)] = exp(—DIFF) x DAIE[C(y0)] X DBIE[C(y)] . (3.10) 


This method of multiplying the Airy and Bairy functions together allows for much larger 
values of C(y) to be used. Thus most of the evanescent waves could be included. 
Underflow errors did not offer the same problems that overflow errors caused since the 
machine handled underflow errors without termination of the computer program 


LSVOCN. However, small inaccuracies from underflow errors had to be expected. 
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These same techniques can be used to solve Eq.(2.76) and Eq.(2.77). Overflow 
errors were not a problem for these two equations since the Bairy function was not part 
of the solution. Graphical presentations of the Airy and Bairy functions are included in 
Figures 7 and 8. Tabulated values of the Airy and Bairy functions, the exponentially 
scaled Airy and Bairy functions, and EDAIE(x) and EDBIE(x) are included in Table 1. 
The functions EDATE(x) and EDBIE(x) are the exponentially scaled Airy and Bairy 
functions multiplied by the appropriate exponential factor that eliminates the exponential 
scaling. These values were also tabulated to ensure that the numerical method 
represented by Eq.(3.9) and Eq.(3.10) would produce the same values as if the 
exponentially scaled functions were not used. This can be observed in Table 1 by 
noticing that DAI(X) = EDAIE(X) and that DBI(X) = EDBIE(X). 

To solve for the velocity potential, Eq.(3.6) must theoretically be integrated from 
zero to infinity. In practice, the upper limit of integration need not be infinity, but rather 
a sufficiently large number so that the integrand becomes negligible in value, in other 
words, all significant evanesent waves are included. A method was developed to 
automatically generate the upper limit of integration. This method will be discussed next. 


We begin with Eq.(2.45), rewritten here for convenience: 
2 
C(y) =—-(01) 3450). (3.11) 
Equation (3.11) can rewritten as 


Cy) =-—_L— Fy - #), : (3.12) 


f(aik3)? 


and upon rearranging terms and solving for k,, we obtain 


kp=+ [e2) + [ fans)? ko ‘ (3.13) 











X _ DAIE(X) 


0.0E+00 0.355028E+00 
0.1E+02 0.158124E+00 
0.1E+03 0.891969E-01 
0.1E+04 0.501642E-01 
O.1E+05 0.282095E-01 
0.1E+06 0.158634E-01 
O.1E+07 0.892062E-02 I a tata aaa alas icicle agai giana 
0.1E+08 0.501643E-02 
0.1E+09 0.282095E-02 
O.1E+10 0.158634E-02 





Airy and Exponentially Scaled Airy Functions 






















ie" > DANG) 
"0" > DAIE() 
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O.1E+16 
O.E+17 
0.1E+18 
O0.1E+19 
0.1E+20 
0.1E+21 
0.1E+22 
0.1E+23 
0.1E+24 
O.1E+25 
0.1E+26 
0.1E+27 
0.1E+28 
0.1E+29 
0.1E+30 
0.1E+31 
0.1E+32 
0.1E+33 
0.1E+34 
0.1E+35 
0.1E+36 
0.1E+37 
0.1E+38 
O.1E+39 
0.1F +40 
O.1E+4) 
0.1E+42 


0.501643E-04 
0.282095E-04 
0.158633E-04 
0.892062E-05 
0.501643E-05 
0.282095E-05 
0.158634E-05 
0.892062E-06 
0.501643E-06 
0.282095E-06 
0.158634E-06 
0.892062E-07 ie Soe 
0.501643E-07 
0.282095E-07 
0.158634E-07 “=e 
0.892062E-08 
0.501643E-08 
0.282095E-08 
0.158634E-08 

0.892062E-09 5) ee 
0.501643E-09 
0.282095E-09 
0.158634E-09 4b----- 
0.892062E-10 
0.501643F-10 
0.282095E-10 
0.158634E-10 





O.1F+11  0.892062K-03 sod ak 2 aa cee ial a ocala aes ae texas 
O.1E+12 0.501643E-03 *! ' ‘ i ' 1 
O.1E+13 0.282095E-03 7 : 
O.1E+14 0.158634E-03 2K : : : 
O.1E+15 0.892062E-04 Ap----- {o----d-----4 jenn nee beeen po----feonnn}e----de---- 











op B------------- 


ge B---- eee eee 












































Pees 
















3 
$ 









1‘ 
) 
5 
t 
t 
t 
1 
1‘ 
‘ 
' 
1 
‘ 
' 
' 
' 
‘ 
4 
‘ 
’ 
4 
1 
! 
! 
' 
1 
---4 
t 
' 
1 
1 
1 
eee eee 












ee eee. Pee 






































ee ne Senay 
ee eee eee Ss 






















Tee atelier en ae 


On Oe ee 
lay eee) eee 


fo sete TUS: 








Rob---- 











Figure 7. Tabulated values of DAIE(x) and plots of DAI(x), DAI(-x), DAIE(x), and EDAIE(x). 
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Bairy and Exponentially Scaled Bairy Functions 

















0.0E+00 0.614927E+00 
0.1E+02 0.318340E+00 
0.1E+03 0.178431E+00 
0.1E+04 0.100329E+00 
0.1E+05 0.564190E-01 
O.1E+06 0.317267E-01 
O0.1E+07 0.178412E-01 
0.1E+08 0.100329E-01 
0.1E+09 0.564189E-02 
O.1E+10 0.317267E-02 
0.1E+11  0.178412E-02 
0.1E+12 0.100329E-02 
0.1E+13 0.564189E-03 
O.1E+14 0.317267E-03 
Q.1E+15 0.178412E-03 
O.1E+16 0.100329E-03 
O.1E+17 0.564190E-04 
O.1E+18 0.317267E-04 
O0.1E+19 0.178412E-04 
0.1E+20 0.100329E-04 
0.1E+21 0.564190E-05 
O.1E+22 0.317267E-05 
0.1E+23 0.178412E-05 
0.1E+24 0.100329E-05 
O.1E+25 0.564190E-06 
0.1E+26 0.317267E-06 
O0.1E+27 0.178412E-06 
0.1E+28 0.100329E-06 
O.1E+29 0.564190E-07 
0.1E+30 0.317267E-07 
0.1E+31 
0.1E+32 0.100329E-07 
0.1E+33 0.564189E-08 
0.1E+34 0.317267E-08 
QO.1E+35 0.178412E-08 
0.1E+36 0.100329E-08 
0.1E+37 0.564189E-09 
0.1E+38 0.317267E-09 
0.1E+39 0.178412E-09 
0.1E+40 0.100329E-09 


0.1E+41 
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Figure 8. Tabulated values of DBIE(x) and plots of DBI(x) 
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DAI(X) DAIE(X) EDATE(X) 
0.355028E+00 0.355028E+00 0.355028E+00 
0.329203E+00 0.336217E+00 0.329203E+00 
0.303703E+00 0.322363E+00 0.303703E+00 
0.278806E+00 0.311084E+00 0.278806E+00 
0.254742E+00 0.301541E+00 0.254742E+00 
0.231694E+00 0.293277E+00 0.231694E+00 
0.209800E+00 0.286000E+00 0.209800E+00 
0.189162E+00 0.279513E+00 0.189162E+00 
0.169846F+00 0.273670E+00 0.169846E+00 
0.151887E+00 0.268364E+00 0.151887E+00 
0.135292E+00 0.263514E+00 0.135292E+00 
0.120049E+00 0.259052E+00 0.120049E+00 
0.106126E+00 0.254928E+00 0.106126E+00 
0.934746E-01 0.251098E+00 0.934746E-01 
0.820380E-01 0.247527E+00 0.820380E-01 
0.717494E-O1 0.244185E+00 0.717495E-01 
0.625369E-01 0.241048E+00 0.625369E-01 
0.543248E-01 0.238094E+00 0.543248E-01 
0.470362E-0!1 0.235306E+00 0.470362E-01 
0.405944E-01 0.232667E+00 0.405944E-01 
0.349241E-01 0.230165E+00 0.349241E-01 
0.299526E-01 0.227786E+00 0.299526E-01 
0.256104E-01 0.225522E+00 0.256104E-01 
0.218320E-01 0.223361E+00 0.218320E-01 
0.185561E-O1 0.221297E+00 0.185561E-01 
0.157259E-01 0.219322E+00 0.157259E-01 
0.132893E-01 0.217429E+00 0.132893E-0] 
O.11198SE-01 0.215613E+00 0.111985E-01 
0.941050E-02 0.213868E+00 0.941051E-02 
0.78863E-02 0.212189E+00 0.788631E-02 
0.659114E-02 0.210572E+00 0.659114E-02 
0.549399E-02 0.209013E+00 0.549399E-02 
0.456744E-02 0.207509E+00 0.456744E-02 
0.378729E-02 0.206056E+00 0.378729E-02 
0.313234E-02 0.204651E+00 0.313234E-02 
0.258410E-02 0.203292E+00 0.258410E-02 
0.212648E-02 0.201976E+00 0.212648E-02 
0.174557E-02 0.200700E+00 0.174557E-02 
0.142939E-02 0.199462E+00 0.142939E-02 
0.116765E-02 0.198261E+00 0.116765E-02 
0.951564E-03 0.197095E+00 0.951564E-03 
0.773630E-03 0.195961E+00 0.773629E-03 
0.627496E-03 0.194859E+00 0.627496E-03 
0.507787E-03 0.193786E+00 0.507787E-03 
0.409974E-03 0.192741E+00 0.409974E-03 
0.330250E-03 0.191724E+00 0.330250E-03 
0.265432E-03 0.190732E+00 0.265432E-03 
0.212861E-03 0.189765E+00 0.212861E-03 
0.170326F:-03 0.188822E+00 0.170325E-03 
0.135992E-03 0.187901E+00 0.135992E-03 
0.108344E-03 0.187002E+00 0.108344E-03 
0.861324E-04 0.186124E+00 0.861324E-04 
0.683285E-04 0.185265E+00 0.683285E-04 
0.540905E-04 0.184426E+00 0.540905E-04 
0.427299E-04 0.183605E+00 0.427298E-04 
0.336853E-04 0.182802E+00 0.336853E-04 
0.265006E-04 0.182015E+00 0.265006E-04 
0.208058E-04 0.181246E+00 0.208058E-04 
0.163017E-04 0.180491F+00 0.16301 7E-04 
0.127471E-04 0.179753E+00 0.127471E-04 
0.994769E-05 0.179028E+00 0.994769E-05 
0.774773E-05 0.178318E+00 0.774773E-05 
0.602246E-05 —0.177622E+00 0.602246E-05 
0.467226E-05 0.176939E+00 0.467226E-05 


DBI(X) DBIE(X) EDBIE(X) 
0.00 0.614927E+00 0.614927E+00 0.614927E+00 
0.10 0.659862E+00 0.646096E+00 0.659862E+00 
020 0.705464E+00 0.664628E+00 0.705464E+00 
0.30 0.752486E+00 0.674409E+00 0.752486E+00 
0.40 0.801773E+00 0.677338E+00 0.801773E+00 
0.50 0.854277E+00 0.674892E+00 0.854277E+00 
0.60 0.911063E+00 0.668324E+00 0.911063E+00 
0.70 0.973329E+00 0.658708E+00 0.973329E+00 
0.80 0.104242E+01 0.646954E+00 0.104242E+01 
0.90 O.11I987E+01 0.633817E+00 0.111987E+01 
1.00 0.120742E+01 0.619912E+00 0.120742E+01 
1.10 0,130707E+01 0.605721E+00 0.130707E+01 
120 0.142113E+01 0.591614E+00 0.142113E+01 
1.30 0.155228E+01 0.577859E+00 0,155228E+0) 
1.40 0.170366E+01 0.564646E+00 0.170366E+01 
1.50 0.187894E+01 0.552094E+00 0.187894E+01 
1.60 0.208247E+01 0.540272E+00 0.208247E+01 
1.70 0.231941E+01 0.529208E+00 0.231941E+01 
1.80 0.259587E+01 0.518898E+00 0.259587E+01 
1.90 0.291918E+01 0.509320E+00 0.291918E+01 
2.00 0.329809E+01 0.500437E+00 0.329809E+01 
2.10 0.374315E+01 0.492203E+00 0.374315E+01 
2.20 0.426704E+01 0.484567E+00 0.426704E+01 
2.30 0.488506E+01 0.477480E+00 0.488506E+01 
2.40 0.561577E+01 0.470890E+00 0.561577E+01 
2.50 0.648166E+01 0.464750E+00 0.648166E+01 
2.60 0.751009E+01 0.459017E+00 0.751009E+01 
2.70 0.873439E+01 0.453648E+H00 0.873439E+01 
2.80 0.101953E+02 0.448608E+00 0.101953E+02 
2.90 0.119426E+02 0.443863EH0 0.119426E+02 
3.00 0.140373E+02 0.439384E+00 0.140373E+02 
3.10 0.165547E+02 0.435146E+00 0.165547E+02 
320 0.195870E+02 0.431125E+00 0.195870E+02 
3.30 0.232483E+02 0.427301E+00 0.232483E+02 
3.40 0.276796E+02 0.423657E+00 0.276796E+02 
3.50 0.330555E+02 0.420177E+00 0.330555E+02 
3.60 0.395927E+02 0.416848E+00 0.395927E+02 
3.70 0.475607E+02 0.413656E+00 0.475607E+02 
3.80 0.572954E+02 0.410592E+00 0.572954E+02 
3.90 0.692160E+02 0.407646E+00 0.692160E+02 
4.00 0.838471E+02 0.404809E+00 0.838471E+02 
4.10 0.101846E+03 0.402075E+00 0.101846E+03 
420 0.124038E+03 0.399435E+00 = 0.124038E+03 
4.30 0.151462E+03 0.396884E+00 0.151462E+03 
4.40 0.185428E+03 0.394417E+00 0.185428E+03 
4.50 0227588E+03 0.392027E+00 0.227588E+03 
4.60 0.280036E+03 0.389712E+00 0.280036E+03 
4.70 0.345426E+03 0.387466E+00 0.345426E+03 
4.80 0.427126E+03 0.385286E+00 0.427126E+03 
4.90 0.529425E+03 0.383168E+00 0.529425E+03 
5.00 0.657792E+03 0.381109E+00 0.657792E+03 
5.10 0.819209E+03 0.379105E+00 0.819210E+03 
5.20 0.102261E+04 0.377155E+00 0.102262E+04 
5.30 0.127947E+04 0.375256E+00 0.127947E+04 
5.40 0.160448E+04 0.373405E+00 0.160448E HM 
5.50 0201658E+04 0.371600E+00 0.201658E+04 
5.60 0.254018E+04 0.369839E+00 0.254018E+04 
5.70 0.320680E+04 0.368120E+00 0.320680E+04 
5.80 0,405720E+04 0.366441E+00 0.405720E+04 
5.90 0.514421E+04 0.364800E100 0.514422E +04 
6.00 0.653645E+04 0.363197E+00 0.653645E+04 
6.10 0.832309E+04 0.361629EH00 0.832309E +04 
620 0.106204E+05 0.360095E+00 0.106204E+05 
630 0.135799E+05 0.358593E+00 0.135800E+05 

















































































































Table 1. Tabulated values for DAI(x), DAYE(x), EDATE(x), DBI(x), DBIE(x), and EDBIE(x). 
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Before the exponentially scaled Airy and Bairy functions were used, the largest 
positive value of the argument for the Bairy function that could be evaluated by the 
computer without causing overflow errors was C(y) = C.,,, = 40. By using the 
exponentially scaled Airy and Bairy functions, it was discovered that C__, could take on 
extremely large values, for example, C_,, = 10” - 10’°. Substituting C__. into Eq.(3.13) 


for C(y) resulted in the following equation for the upper limit of integration. 


er leq) ef [ aaa’)? = (3.14) 


Therefore, the upper limit of integration used in Eq.(3.6) could automatically be adjusted 
to the maximum possible value regardless of the frequency selected. However, due to the 
oscillatory nature of the Bessel function and the dependence of its argument on k, and r, 
the integrand of Eq.(3.6) becomes increasingly difficult to integrate numerically as k, 
increases (for a fixed r), because as the wave number of the Bessel function increases, the 
function becomes more highly oscillatory. To ensure accurate results from the computer 
code, the region of integration had to be divided into subintervals in performing the 
numerical integration. This method enabled LSVOCN to calculate increasingly many 
evanescent waves providing more accurate results. 
B. COMPARISON OF LSVOCN WITH THE TIME-HARMONIC 
FREE-SPACE GREEN'S FUNCTION FOR A HOMOGENOUS OCEAN 
MEDIUM 

1. Green's Function Test Case One - Vary the Source Sound-Speed 

The first test case involved varying the source, operating at frequencies of 
1000 Hz, 250 Hz, and 50 Hz, sound-speed while maintaining the surface sound-speed 
fixed, using both positive and negative gradients for the square of the index of refraction. 
Each case had an initial source sound-speed of 1515 m/s for the negative gradient of the 
square of the index of refraction, and 1485 m/s for the positive gradient case. The source 


and receiver depths were held constant at 150 meters, with the receiver located one meter 








in horizontal range from the source. For each successive run of the test case, the source 
sound-speed was decreased/increased by a factor of 3 m/s until 1501 m/s and 1499 m/s 
was reached, respectively, while maintaining the surface sound-speed constant at 1500 
m/s. Thus, the source sound-speed was never allowed to take on a value of 1500 m/s, 
ensuring that the Airy function solution of LSVOCN, which assumed a sound-speed 
profile where the square of the index of refraction was a linear function of depth, would 
remain valid while at the same time forcing the modeled ocean medium to become more 
homogeneous. If the derived Airy function solution and the corresponding computer code 
in LSVOCN were correct, then as the source sound-speed approached the surface 
sound-speed, the velocity potential (magnitude and phase) values calculated by LSVOCN 
would converge approximately to the velocity potential (magnitude and phase) values 
predicted by Eq.(3.4) and (3.5), respectively. Because we were making calculations at the 
source depth, the DIFF term given by Eq.(3.9) was equal to zero, which made the 
exp(-DIFF) term in Eq.(3.10) equal to one. This condition caused the oscillations of the 
integrand in Eq.(3.6) to decay very slowly, thus requiring longer integration causing an 
increase in CPU time. A summary of these test cases is included as Tables 2 through 7. 

The Delta's and percentage differences included in Tables 2 through 7 were 
calculated as follows: 


1. Delta magnitude = LSVOCN magnitude - Green's function magnitude 
2. Delta phase = LSVOCN phase - Green's function phase 


. ° LSVOCN mag - Green’s function mag 
0 eS 0 
3. % Difference magnitude = ——————_____+* ; x 100% 


LSVOCN phase - Green’s function phase 
Green’s function phase 


4. % Difference phase = x 100% 

The data in Tables 2 through 7, show that as the source sound-speed approaches 
the surface sound-speed, making the modeled ocean medium more homogeneous, the 
calculated velocity potential (magnitude and phase) values from LSVOCN approach the 
velocity potential (magnitude and phase) values predicted by the Green's function - our 


first validation. 
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SOURCE FREQUENCY = 1000Hz 
SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s 


Phase Magnitude Phase Delta Delta % Difference 
Green's Func] Green's Func | LSVOCN LSVOCN Magnitude Phase Magnitude 
Deg. m/s Deg. m/s Deg. 


[1512 0.0795775_] 301.90476 
| 1509 0.0795775 J 301.43141  0.0788504 [301.559 -0.0007271 | 0.12759_ | -0.9137005 | 0.042328 _| 
| 1506 J 0.0795775 J 300.95618 ] 00789223 J 301.039 | -0.0006552 0.08282 J -0.8233483 | 0.027519 _| 


| 1501 | 00795775 J 300.15989 | 0.0791721 0.01189 


Table 2. Comparison of Green's Function with LSVOCN for Test Case 1 
(Vary the Source Sound-Speed) - n?(y) with Negative Gradient. 





SOURCE FREQUENCY = 250Hz 
SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s 


Sound Speed§ Magnitude Phase Magnitude Phase Delta Delta % Difference | % Difference 
Source Green's Func | Green's Func |. LSVOCN LSVOCN Magnitude Phase Magnitude Phase 
m/s m/s Deg. mss Deg. m/s Deg. 


S15 0.0795775_ J 120.59406 J 0.0791376 | 120.992 fl -0.0004399 T_ 0.39794 -0.5527944 | 0.3299831_ | 
| 1512 | 0.0795775 J 120.47619 | 0.0791609 | __ 120.83 f{ -0.0004166 | _0.35381_ | -0.5235148 | 0.2936763_| 
| 1509 0.079575 | 120.35785 | 0.0792156_] 120.704 -0.0003619 [0.34615 _ | -0.4547768 }0.2876007_ | 
| 1506 J 0.0795775_f_120.23904 | 0.07923 J 120.502_ | -0.0003475 J 0.26296 _ f -0.4366812 J 0.2186977_| 
|__ 1503 0.079575 | 120.11976 | 0.079377 J 120.422] -0.0001998 J 0.30224 -0.251076 I 0.2516156_| 
| 1501 0.0795775 I 120.03997 J 0.0794475 T120.263 T -0.00013_ F0.22303 ft -0.1633628 J 01857964 


Table 3. Comparison of Green's Function with LSVOCN for Test Case 1 
(Vary the Source Sound-Speed) - n?(y) with Negative Gradient. 





SOURCE FREQUENCY = 50Hz 
SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s 


Sound Speed J Magnitude Phase Magnitude Phase Delta Delta % Difference 
Source Green's Func j Green's Func LSVOCN LSVOCN Magnitude Phase Phase 
ms mis * Deg. nr/s Deg. m/s Deg. 


0.24719 
0.0001308 
0.0002397 


Table 4. Comparison of Green's Function with LSVOCN for Test Case 1 
(Vary the Source Sound-Speed) - n?(y) with Negative Gradient. 

















SOURCE FREQUENCY = 1000Hz 
SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s 


Phase Magnitude Phase Delta Delta % Difference | % Difference 
Green's Func { LSVOCN LSVOCN Magnitude Phase Magnitude Phase 
Deg. m’/s Deg. m/s Deg. 


















Sound Speed 
Source 
m/s 


Magnitude 
Green's Func 
m/s 





0 
1491 f0.0795775 J 29855131 | o07es4es | 298.647 -0.0007291 | 0.09569 | 09162138 | 0.030510 1 
[1494] 0.0795775 J 299.03614 | 0.0789214 | 299.099 | -0.0006561 | 0.06286 J -0.8244793 
11499} 0.0795775 ]_ 29983989 | 0.079143 0.00089 









Table 5. Comparison of Green's Function with LSVOCN for Test Case 1 
(Vary the Source Sound-Speed) - n?(y) with Positive Gradient. 






SOURCE FREQUENCY = 250Hz 
SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s 


Sound Speed] Magnitude Phase Magnitude Phase Delta Delta % Difference | % Difference 
Source Green's Fanc } Green's Func | LSVOCN LSVOCN Magnitude Phase Magnitude Phase 
nas m/s Dep. ws Deg. m/s Deg. 


1485 0.0795775_ 4 119.39394 J 0.0791194 J 119.774 J -0.0004581 | 0.38006 {f -0.5756652 | 0.3183244 
1488 0.0795775 J 119.51613 | 0.0791708 | 119.885 0.36887 | -0.5110741 | 0.3086362 
0795775 J _119.63783 J 0.0791928 F_ 119.958 f -0.0003847 | 0.32017 | -0.4834281 | 0267616 

0 





nN 


|_i4oi fo 1 
| 1494 0.0795775 T119.75904 0.2513046 
se en 0795775 J 119.95997 F 0.0794502 [120.185] -0.0001273 J 0.22503 J -0.150960R ] O.1875876 | 


Table 6. Comparison of Green's Function with LSVOCN for Test Case 1 
(Vary the Source Sound-Speed) - n?(y) with Positive Gradient. 


SOURCE FREQUENCY = 50Hz 
SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s 


Magnitude Delta % Difference 
LSVOCN Magnitude Phase 
m’/s m/s 


0.1132408 


0.0946862 
0.0660805 


Table 7. Comparison of Green's Function with LSVOCN for Test Case | 
(Vary the Source Sound-Speed) - 7?(y) with Positive Gradient. 








2. Green's Function Test Case Two - Vary the Horizontal Range 

The second test case involved varying the horizontal range while maintaining the 
source and receiver depths fixed at 150 meters. For this case the surface and source 
sound-speeds are not varied. If the receiver is fixed at the source depth (150 meters) and 
the horizontal range is increased, then the magnitude of the velocity potential should 
decrease by a factor of 1/4nR, independent of frequency, as predicted by Eq.(3.4). 
Furthermore, it should be observed that the phase varies as predicted by Eq.(3.5): 


£oy=4( 21 - 2. (3.15) 


The above observations are reasonable assumptions to make even though the modeled 
ocean medium is not homogeneous. As long as the receiver is maintained at the same 
depth as the source, then the sound speed is constant, which simulates a homogeneous 
ocean medium. 

To conduct the second Green's function test case, the receiver depth and source 
depth are both fixed at 150 meters. The surface sound-speed is fixed at 1500 m/s and the 
source sound-speed is fixed at 1515 m/s for a negative gradient for the square of the 
index of refraction, and 1485 m/s for a positive gradient for the square of the index of 
refraction. Horizontal range values selected were 2.0, 3.0, and 5.0 meters. These values 
were selected based on the amount of CPU time required to calculate the velocity 
potential (magnitude and phase) for each selected receiver depth, horizontal range pair. 
Because we were making calculations at the source depth, the D/FF term given by 
Eq.(3.9) was equal to zero, which made the exp(-D/FF) term in Eq.(3.10) equal to one. 
This condition caused the oscillations of the integrand in Eq.(3.6) to decay very slowly, 
thus making integration more difficult resulting in more CPU time being required. In 
addition, as k increases (for a fixed r) the zeroth-order Bessel function in Eq.(3.6) 
oscillates more rapidly, making the numerical integration more difficult. Any calculation 
requiring more than 60 minutes of CPU time was stopped by the computer. The test case 


results are summarized in Tables 8 and 9. Calculations that required greater than 60 
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minutes of CPU time are denoted by: >60 min. The delta's and percentage differences 
were calculated in the same manner as discussed previously. The data in Tables 8 and 9 
show that as the horizontal range increases while maintaining the receiver at the source 
depth, the velocity potential (magnitude and phase) calculated by LSVOCN behaves 
according to the Green's function for a free-space homogeneous ocean medium as 


predicted by Eq.(3.4) and Eq.(3.5) with S, = 1 - our second validation. 


SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s SOURCE SOUND-SPEED = 1515 m/s 





0: 
0.01591559120.5940690.0157484 121.152 | -0.0001671 0.55794] -1.0499199 0.4626596 





Table 8. Comparison of Green's Function with LSVOCN for Test Case 2 
(Vary Horizontal Range ) - n?(y) with Negative Gradient. 












SOURCE DEPTH = RECEIVER DEPTH = 150 METERS 
SURFACE SOUND-SPEED = 1500 m/s SOURCE SOUND-SPEED = 1485 m/s 


Frequency Range LSVOCN Dehta Delta % % 
Hertz meters Phase Magnitude Phase Difference | Difference 
Deg. 6 Deg Magnitude Phase 
/ A 
a 
/ 





on 
N/A 





iA / 

‘A 

ee Te 
a Cd ET eT 

Goo1sea] 02917989 
| 


119.30394 


Table 9. Comparison of Green's Function with LSVOCN for Test Case 2 
(Vary Horizontal Range ) - 7?(y) with Positive Gradient. 
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As a final note on the second Green's function test case, as can be observed from 
Tables 8 and 9, the delta magnitudes all agree out to three significant digits, and in some 
cases better. This fact was consistent with the magnitude differences as calculated in the 
first Green's function test case. Also, for small velocity potential phases, small phase 
errors result in large percentage differences, as expected. 

3. Green's Function Test Case Three - Vary the Receiver Depth 

The third Green's function test case involved varying the receiver depth above 
and below the source depth, while maintaining the source and surface sound-speeds fixed 
at 1515/1485 m/s and 1500 m/s, respectively. As was stated before, the modeled ocean 
medium is not homogeneous; however, for the previous two Green's function tests, the 
receiver depth was maintained at the source depth and the source sound-speed was varied 
to more closely match the surface sound-speed (Test Case 1) or maintained constant 
(Test Case 2) which simulated a homogeneous ocean medium. For this case, the receiver 
was moved slightly above and below the source depth, and as a result, the Green's 
function results are less valid. As can be seen from Figure 9, even though the receiver is 
not located at the source depth, the actual change in the sound speed at the receiver is 
very small. Because the ocean medium acts very much like a homogeneous medium 
when receivers and sources are "close", the velocity potential (magnitude and phase) 
computed using the free-space Green's function for a homogeneous ocean medium is a 
valid approximation. This test case allowed us to validate the remaining part of the 
computer code in LSVOCN for a source and receiver at different depths, prior to 
comparing it with the RRA Algorithm. This test case was conducted with a source 
operating at 1000 Hz (excluding cases for the source operating at 250 Hz and 50 Hz since 
validation at one frequency was adequate). The third Green's function test case is 


summarized in Table 10. The test case was conducted four times: 


1. For the square of the index of refraction with a negative gradient and a 
receiver depth of 150.5 meters and a horizontal range of 0.8660254 meters. 


2. For the square of the index of refraction with a negative gradient and a receiver 
depth of 149.5 meters and a horizontal range of 0.8660254 meters. 








3. For the square of the index of refraction with a positive gradient and a receiver 
depth of 150.5 meters and a horizontal range of 0.8660254 meters. 


4. For the square of the index of refraction with a positive gradient and a receiver 


depth of 149.5 meters and a horizontal range of 0.8660254 meters, 


Region I: Receiver above the source 


Receiver Depth 


R=1 met 
149.5 meters said 


r= 0.8660254 meters 


Source depth 
150 meters 


Receiver Depth 
150.5 meters 


Region II: Receiver below the source 





Figure 9. Geometry for Green's Function Test Case 3 (Vary Receiver Depth). 


The surface sound-speed is maintained at 1500 m/s. The source sound-speed is 
maintained at 1515 m/s when the square of the index of refraction has a negative 
gradient, and 1485 m/s when the square of the index of refraction has a positive gradient. 
Calculations are based on Eq.(3.4) and Eq.(3.5) where R and c are treated as variables. 
As can be observed from Table 10, the change in the sound speed at the location of the 
receiver is very small, making our assumptions for this test case reasonable. Analysis of 
the data presented in Table 10 show that LSVOCN calculated values for the velocity 
potential (magnitude and phase) are approximately equal to those predicted by the 
Green's function - our third validation. The Delta's and percentage differences were 


calculated in the same manner as discussed before. 
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SURFACE SOUND-SPEED = 1500 M/S 
r= 0.8660254 METERS, FREQ = 1000 Hz 


Source | Receiver Sound | Receiver} Green's Green's | LSVOCN J LSVOCN Delta Delta % Diff % Diff 
Sound Speed Depth Func Fune Mag Phase Mag Phase Mag Phase 
Speed m/s meters | Magnitude | Phase m/s deg. m/s deg. 

m/s m/s deg. 


| 1515 | 1514.9492501 [149.5 | 
1485 | 1485.0492549 
Table 10. Comparison of Green's Function with LSVOCN for Test Case 3 
(Vary the Receiver Depth) - n?(y) with Negative and Positive Gradients. 







C. COMPARISON OF LSVOCN WITH THE RECURSIVE RAY ACOUSTICS 
ALGORITHM 


1. Test Preparation 

The RRA Algorithm calculates the magnitude in Pascals and phase in degrees of 
the acoustic pressure along a ray path. The source is modeled as a time-harmonic, 
omnidirectional point source, with its strength represented by a source level (SL) value in 
decibels relative to one micro-pascal (rms). The "Airy function solution" is in terms of 
velocity potentials and is based on a time-harmonic, omnidirectional point source 
modeled by a unit amplitude impulse function. Thus, before the magnitudes and phases 
from the RRA Algorithm and LSVOCN (the "Airy function solution") could be 
compared, the velocity potential expressions were transformed into acoustic pressure 


expressions using the following equation for time-harmonic fields [Ref. 2, Chapter 2]: 


PAT, Y) = —j20fpoo/r,y) (3.16) 


where f is the frequency of the source in Hz, and Po is the ambient or equilibrium 
density of the fluid medium in kg/m’. 

In addition, the magnitudes from LSVOCN had to be made equivalent to the 
magnitudes computed by the RRA Algorithm by relating the source level in decibels to 
the amplitude of the impulse function used to model the omnidirectional point source in 


LSVOCN. This was accomplished by deriving a magnitude scale factor that relates a 








source level in decibels to the source strength of a spherical sound source operating in the 


monopole mode of vibration as the radius of the sphere approaches zero [see Eq.(3.3)]. 


Thus, substituting Eq.(3.3) into Eq.(3.16) yields 


KR 
PAR , 0, W) = +i f poSo ,R>0. 
The magnitude of Eq.(3.17) is given by 


lpAR , 0, w)| = es, /R>O. 


Thus, the peak acoustic pressure, P, , at R equal to one meter is 


Po= lpAR 9, Wr = es, : 


Therefore, the source strength, S, , in m’/s, is given by 


(3.17) 


(3.18) 


(3.19) 


(3.20) 


The peak acoustic pressure, P,, in Pascals [Ref. 2, Chapter 1] due to a given source level 


in decibels is given by: 


Po = (2 P02) 


(3.21) 


where P,_-is the reference pressure equal to one micro-pascal (rms), and SL is the source 


level of the source in decibels relative to P_,. Since the source level used by the RRA 


Algorithm is a known constant (180.0 dB), the peak acoustic pressure can be calculated 


using Eq.(3.21). This result can then be substituted into Eq.(3.20) allowing the source 
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strength to be computed. The resulting source strength value, S,, is the magnitude scale 
factor that was used to multiply the output from LSVOCN. 

In summary, the conversion process to make the LSVOCN magnitude and phase 
equivalent to that computed by the RRA Algorithm is as follows [multiply the right-hand 
side of Eq.(3.16) by Eq.(3.20)]: 


pAr.)) = Pont poe "Ar,y) = H4RPoOy ry) (3.22) 


where P, is given by Eq.(3.21). Equation (3.22) was added as FORTRAN computer code 
to LSVOCN so that the conversion process would be taken care of automatically. 

2. Comparison of the RRA Algorithm with LSVOCN when n’(y) has a 

Negative Gradient 

The RRA Algorithm was used to simulate various ray paths propagated through 
the ocean medium when the square of the index of refraction was a linear function of 
depth with a negative gradient, corresponding to a positive gradient sound-speed profile. 
Two ray paths identified by the initial launch angle, B, , were selected for each of the test 
frequencies, 1000 Hz, 250 Hz, and 50 Hz. The sound-speed profile and the ray 
propagation path for each ray is illustrated in Figure 10. Launch angles were selected so 
that the ray paths did not interact with the ocean surface or bottom, thus ensuring that the 
magnitude and phase calculated by the RRA Algorithm were for a free-space propagation 
problem. Test points along each ray path were selected. The depth (receiver depth) and 
the horizontal range (HRNG) corresponding to these points were recorded and used as 
inputs to LSVOCN. The comparison between the RRA Algorithm and LSVOCN 
magnitude and phase values are summarized in Tables 11 through 13, where the asterisk 
"*" means that the entire row corresponds to data at the location of a turning point. The 
calculated values in Tables 11 through 13 were obtained as follows: 


1. Delta magnitude = RRA Magnitude - LSVOCN Magnitude 


2. Delta phase = RRA Phase - LSVOCN Phase 








Ww 


. Initial phase offset = Delta phase of the first point for each ray 


4. Adj RRA Phase = (RRA Phase - Initial phase offset) 
a) Add 360°if Adj RRA Phase angle is negative 
b) Subtract 360° if Adj RRA Phase angle is greater than 360° 
c) Do not add or subtract if Adj RRA Phase angle is positive and less than 360° 


RRA Magnitude - LSVOCN Magnitude 
LSVOCN Magnitude 


ws 


. % Difference Magnitude = x 100% 


Adj RRA Phase - LSVOCN Phase 
LSVOCN Phase 


an 


. % Difference Phase = x 100% 


Analysis of the data if Tables 11 through 13 indicate that the magnitudes and 
phases calculated by the RRA Algorithm match reasonably well with the magnitudes and 
phases calculated by LSVOCN. Note that there is an initial phase difference that exists 
between the RRA Algorithm and LSVOCN. The RRA Algorithm begins with zero phase 
and calculates the phase thereafter based on travel time as a ray propagates. In order to 
correct for this initial phase difference, each RRA phase angle after the first data point on 
each ray is adjusted by an amount equal to the delta phase of the first data point. After 
the RRA phase has been adjusted, a valid comparison between the RRA phase and 
LSVOCN phase can be performed. It can further be observed that the percentage phase 
error tends to increase as the ray propagates. This seems reasonable since small errors 
made in calculating travel time and, hence, phase in the beginning accumulate as the ray 
propagates further from the source. These phase errors could possibly be reduced by 


selecting a smaller arc length step size used by the RRA Algorithm. 
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SURFACE SOUND-SPEED = 1500 m/s, SOURCE-SOUND SPEED = 1515 m/s 
SOURCE DEPTH = 150 m, n’(y) WITH NEGATIVE GRADIENT 
FREQUENCY = 1000 HZ, "*" TURNING POINT DATA 


RRA | RRA |LSVOCN|LSVOCN | DELTA | DELTA] Adj RRA | %DIFF | %DITr 
Mag Phase Mag Phase Mag Phase Phase Mag Phase 
Pascals | Degree | Pascals | Degrees | Pascals Degrees | Degrees 

1414.21 | 122.377] 1414.06 | 212.411 | 0.15 0.01061 N/A 


500.178 2,79558 -0.85388 
1.40685 303.34 0.62504 
0.94006 0.76272 
| * 294.58 _[ 2056.996 | 0.68751 1.32725 
241.2 -2.27488 
75.886 3,22457 
|_1__ [498.762 [ 0.28757 [238.727] 0.30947 | 337.044 | -0.0219 | -98.317 | 326.865 | -7.07661 | -3.02008 
piso [1 [aia Tizzs76] aio | 212.154 | 3a [20778 | Na ores] NAT 























116.533, [999.753 | 1.41456 | 39.786 | i.e4o1s | 131.49 |-0.02555] -01.704 | 129.568 | 77417 | ni aoa7s | 
209.696 
Hae 25 .2000.033 | 0.7071 [341.127] 0.66581 [69.119 | 0.04129 | 272.008 | 70.905 | 20147 | 258308 


Table 11. Comparison of RRA Algorithm with LSVOCN Magnitude and Phase Values for different 
Receiver Depth and Horizontal Range Points Located on Two Separate Sound Rays, Identified by Bo. 












SURFACE SOUND-SPEED = 1500 m/s, SOURCE SOUND-SPEED = 1515 m/s 
SOURCE DEPTH = 150 m, n°(y) WITH NEGATIVE GRADIENT 
FREQUENCY = 250 HZ, "*" TURNING POINT DATA 


RRA RRA LSVOCN | LSVOCN | DELTA | DELTA Adj RRA | % DIFF | % DIFF 
Mag Phase Mag Phase Mag Phase Phase Mag Phase 
Pascals Degree Pascals Degrees | Pascals | Degrees Degrees 
b 150.139 | 0.99 [141421 T 300.594 | 1415.07 06077 |N/A 


500.178 1.10857 
| 256.408 1.37561 | 325.296 0.4439 
1500.38 0.90205 0.54563 
* 294.58 |2056.996] 0.68751 | 234.732 | 0.65279 0.90565 
-0.00954 ~4.96258 
f 90.742 [4499.953] 0.31804 352.622 | 0.04319 357,127 1.27757 


|i [498.762] 0.28757 | 59.682 | 0.2389 | 151.482 | 0 04R67 149,734 | 20.37254{ -1.15393 
Be erat ale ates ee ee 








i Receiver | Receiver 





‘ 
oO 


















Table 12. Comparison of RRA Algorithm with LSVOCN Magnitude and Phase Values for different 
Receiver Depth and Horizontal Range Points Located on Two Separate Sound Rays, Identified by Bo. 








SOURCE DEPTH = 150 m, n’(y) WITH NEGATIVE GRADIENT 
FREQUENCY = 50 HZ, "*" TURNING POINT DATA 


Receiver | Receiver | RRA RRA |LSVOCN]| LSVOCN | DELTA | DELTA | Adj RRA % DIFF | % DIFF 
f Depth HRNG Mag Phase Mag Phase Mag Phase Phase Mag Phase 
meters meters Pascals | Degree | Pascals | Degrees | Pascals | Degrees Degrees 
Ee a eee epee ee ee ie, 
1 150.139 | 0.99 | 1414.21 [348.119] 1396.42 270.155 1.27397 


f 211.757 | 500.178 | 2.80646 144.851] 2.7779 | 235.583 | 0.02856 | -90.732 | __234.696__| 1.02811 | -0.37651_| 


-1,59797 
| 116.533 | 999.753 | 1.41456 |343.989| 1.4166 | 75.872 _|-0.00204| 268.117] 74.081 __| -0.14401 
-1.47554 
f 16.029 | 2000.033 | 0.7071_]215.056] 0.67256 | 302.346 | 0.03454 | -87.29 | 305.148 | 5.1356 | 0.92675 F 


Table 13. Comparison of RRA Algorithm with LSVOCN Magnitude and Phase Values for different 
Receiver Depth and Horizontal Range Points Located on Two Separate Sound Rays, Identified by Bo . 





3. Comparison of the RRA Algorithm with LSVOCN when n’(y) has a 
Positive Gradient 

The RRA Algorithm was used to simulate various ray paths propagated through 
the ocean medium when the square of the index of refraction was a linear function of 
depth with a positive gradient, corresponding to a negative gradient sound-speed profile. 
Two ray paths identified by the initial launch angle, B, , were selected for each of the test 
frequencies, 1000 Hz, 250 Hz, and 50 Hz. The sound-speed profile and the ray 
propagation path for each ray is illustrated in Figure 11. Launch angles were selected so 
that the ray paths did not interact with the ocean surface or bottom, thus ensuring that the 
magnitude and phase calculated by the RRA Algorithm were for a free-space propagation 
problem. Test points along each ray path were selected. The depth (receiver depth) and 
the horizontal range (HRNG) corresponding to these points were recorded and used as 
inputs to LSVOCN. The comparison between the RRA Algorithm and LSVOCN 
magnitude and phase values are summarized in Tables 14 through 16, where the asterisk 


"*" means that the entire row corresponds to data at the location of a turning point. The 
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calculated values in Tables 14 through 16 were obtained in the same manner as was done 
for Tables 11 through 13. Analysis of the data in Tables 14 through 16 indicate that the 
magnitudes and phases calculated by the RRA Algorithm match reasonably well with the 
magnitudes and phases calculated by LSVOCN. Note that there is an initial phase 
difference that exists between the RRA Algorithm and LSVOCN. The RRA Algorithm 
begins with zero phase and calculates the phase thereafter based on travel time as a ray 
propagates. In order to correct for this initial phase difference, each RRA phase angle 
after the first data point on each ray is adjusted by an amount equal to the first delta 
phase of the first data point. After the RRA phase has been adjusted, a valid comparison 
between the RRA phase and LSVOCN phase can be performed. It can further be 
observed that the percentage phase error tends to increase as the ray propagates. This 
seems reasonable since small errors in the beginning accumulate as the ray propagates 
further from the source. These phase errors could possibly be reduced by selecting a 
smaller arc length step size used by the RRA Algorithm. Also, for small acoustic 


pressure phases, small errors result in large percentage differences, as expected. 


SURFACE SOUND-SPEED = 1500 m/s, SOURCE SOUND-SPEED = 1485 m/s 
SOURCE DEPTH = 150 m, n*(y) WITH POSITIVE GRADIENT 
QUENCY = 1000 HZ, "*" TURNING POINT DATA 


FRE 
Receiver | Receiver RRA RRA |LSVOCN|LSVOCN| DELTA | DELTA | RRA Adj | %DIFF | % DIFF 
HRNG Mag Phase Mag Phase Mag Phase Phase Mag Phase 
meters Pascals Degree Pascals_ | Degrees {| Pascals | Degrees | Degrees 


207.61 
09.598 -89.189] 210.442] -0.79219 
-88.83 


-0.03631] _-88.168{ 300.866] -3.72002 
7] 5.24143 
0.3858 0.02111 0.693] 5.47391] -90.0259 
0.3614 0 132.296] -12.0228 
0 0. 7 


2397 19.9491 
[207.709 | 1496] 90.133 
2 





Table 14. Comparison of RRA Algorithm with LSVOCN Magnitude and Phase Values for different 
Receiver Depth and Horizontal Range Points Located on Two Separate Sound Rays, Identified by Bo. 
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SURFACE SOUND-SPEED = 1500 m/s, SOURCE SOUND-SPEED = 1485 m/s 
SOURCE DEPTH = 150 m, n’*(y) WITH POSITIVE GRADIENT 
FREQUENCY = 250 HZ, "*" TURNING POINT DATA 


Receiver | Receiver RRA RRA |LSVOCN|/LSVOCN]| DELTA | DELTA |] RRA Adj | % DIFF % DIFF 
Depth HRNG Mag Phase Mag Phase Mag Phase Phase Mag Phase 
meters meters Pascals Degree Pascals |] Degrees | Pascals | Degrees | Degrees 


1149.861{ 0.99 | 1414.21 | 299.394 | 1415.07 | 29.447 | -0.86_| 269.947 | N/A__| -0.06077 | N/A__| 

[202.264 | 4499.27 | 0.31802 | 10.566 | 0.30859 | 91.865 _| 0.00943 | -81.299 | 100.619 | 3.05583 | 9.5292 | 

0.04464 
B=90° 

| 150_[ 1 | 1414.21 [ 299.304 | 1406.07 | 29.774 | 8.14 | 269.62 | N/A | 0.57892 |N/A __| 

158.268 | 499.534 | 2.83107 | 314.164 | 2.80951 | 44.984 | 0.02156 | 269.18 | 44.544 | 0.76739 | -0.97813 | 

t 282.604 | 1999.778 | 0.70719 | 130.724 | 0.67483 | 223.51} 0.03236 | -92.786 | 221.104 | 4.79528} —1.07646_| 


Table 15. Comparison of RRA Algorithm with LSVOCN Magnitude and Phase Values for different 
Receiver Depth and Horizontal Range Points Located on Two Separate Sound Rays, Identified by B. : 





SURFACE SOUND-SPEED = 1500 m/s, SOURCE SOUND-SPEED = 1485 m/s 
SOURCE DEPTH = 150 m, n°(y) WITH POSITIVE GRADIENT 
QUENCY = 50 HZ, "*" TURNING POINT DATA 


FRE 
Receiver} Receiver RRA RRA {LSVOCNILSVOCN] DELTA [| DELTA | RRA Adj | % DIFF % DIFF 
meters meters Pascals Degree Pascals | Degrees | Pascals | Degrees | Degrees 
ee 
149.861 | 0.99 | 1414.21 | 347.879 | 1395.57 | 77.713 | 18.64 | 270.166 | N/A__| 1.33565 | N/A__| 
B90" | 

| sof 1414.21 [347.879 | 1417.16 | 78.123 | -2.95_ | 269.756 | N/A | 0.20816 | N/A 


Table 16. Comparison of RRA Algorithm with LSVOCN Magnitude and Phase Values for different 
Receiver Depth and Horizontal Range Points Located on Two Separate Sound Rays, Identified by B.. 
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IV. SUMMARY AND RECOMMENDATIONS 

This thesis had two primary goals. First, derive and document the solution to the 
three-dimensional inhomogeneous Helmholtz equation for a free-space propagation 
problem when the square of the index of refraction of the ocean medium is a linear 
function of depth, and the source is an omindirectional point source. Given this type of 
index of refraction, the solution to the inhomogeneous Helmholtz equation has an exact 
solution in terms of Airy functions. This exact solution, in terms of Airy functions, was 
then incorporated into the computer program Linear Space-Variant Ocean (LSVOCN) as 
an additional ocean medium transfer function. The second goal was to further test the 
Recursive Ray Acoustics (RRA) Algorithm by comparing the magnitudes and phases of 
the acoustic sound pressure calculated, at various points, along a ray path generated by 
the RRA Algorithm with corresponding magnitudes and phases calculated by LSVOCN 
along the same ray path. 

In Chapter II, the solution to the aforementioned Helmholtz equation was 
carefully derived and documented. This documentation served to record in detail the 
steps that led to the final Airy function solution. 

The next step in the solution process was to convert the theoretical Airy function 
solution to the FORTRAN computer code that was used by LSVOCN. The main 
problems in the conversion process were the many numerical errors that developed. 
There were two main numerical errors encountered, overflow and underflow errors. 
Overflow errors resulted from the evaluation of the Bairy function with large positive 
argument values. The Airy function did not have this problem since it is well behaved for 
both negative and positive argument values. A graphical representation of both the Airy 
and Bairy functions is shown in Figures 7 and 8, respectively. Overflow errors result in 
the termination of the computer code before it has finished a calculation. Thus no answer 
is ever obtained - a serious problem. To prevent overflow errors, we used the 
exponentially scaled Airy and Bairy functions, also shown in Figures 7 and 8. These two 
functions allowed a significant advance towards a working FORTRAN computer code 


that would solve the "Airy function solution" as given by Eq.(2.75) through Eq.(2.78). 





Using the exponentially scaled functions allowed for very large positive values of C(y), 
the argument of the Airy and Bairy functions, to be used. These large positive values of 
C(v) were the result of propagating waves becoming evanescent waves. For the solution 
to be accurate, it is very important to account for these evanescent waves. The 
exponentially scaled functions did not prevent underflow errors. These types of errors do 
introduce some numerical inaccuracies, but do not result in computer program 
termination. 

The last numerical problem to solve was to establish an upper limit of integration 
for Eq.(2.3). Using the computer, it is impossible to literally integrate to infinity. To 
achieve numerical accuracy, the upper limit of integration had to be large enough to 
account for most of the significant evanescent waves. Thus, a method to automatically 
compute this upper limit of integration, regardless of the source frequency selected, was 
developed. 

Before comparing the RRA Algorithm with LSVOCN, we developed test cases to 
validate the Airy function solution and FORTRAN computer code used by LSVOCN. 
This validation was carried out by comparing the magnitude and phase of the velocity 
potential equal to the free-space Green's function of a homogeneous ocean medium with 
the magnitude and phase of the velocity potential computed by LSVOCN. The 
free-space Green's function was chosen because we could make the modeled ocean 
medium used by LSVOCN "simulate" a homogeneous medium without actually being 
made homogeneous. Therefore, the Green's function provided a simple theoretical 
solution that we could compare the calculated solution of LSVOCN with. All three 
Green's function test cases discussed in detail in Chapter II] gave very good results. Now 
we were finally ready to test calculated values from the Recursive Ray Acoustics (RRA) 
Algorithm with those predicted by LSVOCN. 

The RRA Algorithm versus LSVOCN comparison tests were conducted first, 
when the square of the index of refraction had a negative gradient and second, when the 
square of the index of refraction had a positive gradient. The test results indicated that 


the magnitudes from the RRA Algorithm and LSVOCN were in close agreement. In 
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order to be able to compare the RRA Algorithm with LSVOCN, a magnitude scale factor 
was derived that related the source level in decibels used by the RRA Algorithm to the 
amplitude of the impulse function used to model the omnidirectional point source in 
LSVOCN. The magnitudes of LSVOCN were extremely close to the values predicted by 
the RRA Algorithm. 

The phase comparisons between the RRA Algorithm and LSVOCN were also 
quite good. The initial phase offset between LSVOCN and the RRA Algorithm was 
determined by calculating the phase differences, one meter from the source, between the 
RRA Algorithm and LSVOCN, i.e., the phase difference between the RRA Algorithm 
and LSVOCN of the first data point of each ray path. Once this difference was 
determined, it was subtracted from the other phase values, less the first data point, 
calculated by the RRA Algorithm along the same ray path. It was observed that the 
percentage phase error tended to grow as the receiver was moved further from the source. 
This trend makes intuitive sense because it suggests that small phase errors in the 
beginning are propagated and accumulated as the ray travels further from the source. It 
was also observed that several data points resulted in large percentage differences for the 
phase. These large percentage differences, resulting from the variance between small 
phases could possibly be reduced by making the integration step size a function of the 
wave number rather than a fixed quantity. 

Because ray acoustics is an approximate solution of the wave equation, 
characterized as a "high frequency" approximation, we expected to see better and better 
agreement between the RRA Algorithm and LSVOCN as the source frequency increased 
from 50 Hertz to 1000 Hertz. The data in Tables 11 through 16 did not indicate such a 
trend for all cases. In fact, some percentage errors in magnitude and phase are larger for 
the 1000 Hertz case than the 50 Hertz case. The expected results could possibly be 
observed if the integration step size were made a function of the wave number, as 
discussed earlier. The results from all six test cases comparing the RRA Algorithm with 
LSVOCN indicate that the RRA Algorithm was, in general, accurate regardless of the 


frequency selected. 








As a final note, the CPU time required by the RRA Algorithm to compute all of 
the magnitudes and phases for thousands of points along each ray path was typically on 
the order of 30 seconds. Table 17 gives the approximate CPU times for each LSVOCN 
data point calculated. Examination of Table 17 reveals that the RRA Algorithm is a much 
more efficient method to calculate the acoustic pressure compared to the Airy function 
wave solution. The wave solution not only takes considerably more time, but is also more 


susceptible to numerical inaccuracies. 








CPU TIMES PER DATA POINT CALCULATED BY LSVOCN 


B.= 90° B.= 82° B,= 90° B.= 98° 
0 0 0 0 


Square of the index of refraction with Square of the index of refraction with 
Negative gradient Positive gradient 


receiver @ 1 meter | receiver @ I meter | receiver @ 1 meter | receiver @ | meter 







































~59 min ~50 min ~54 min ~22.5 min 
other receiver points | other receiver points | other receiver points{ other receiver points 
~6 min ~6 min ~6 min ~6 min 





















receiver @ 1 meter | receiver @ 1 meter | receiver @ 1 meter 






~ 22 min ~ 16 min ~ 23 min 
other receiver point Jother receiver point | other receiver point 
~5 min ~5.5 min ~5.3 min 
















receiver @ 1 meter 
~26 min 

other receiver point 

~5 min 


receiver @ 1 meter 
~8.2 min 

other receiver point 

~4 min 


receiver @ 1 meter 
~25 min 

other receiver point Jother receiver point 

~5 min ~5 min 


Table 17. CPU Times Per Data Point Calculated by LSVOCN. 


receiver @ 1 meter 
~8 min 









Recommendations for further research include: 


1. To achieve more accurate calculated acoustic pressure values from the RRA 
Algorithm (magnitude and phase) a smaller arc length step size, used by the 
RRA Algorithm, should be selected. Using a smaller arc length step size may 
decrease small phase errors that are propagated and accumulated as the ray 
travels further from the source. 


2. Propagate sound rays, using the RRA Algorithm, such that focal points are 
generated. Compare the RRA Algorithm and LSVOCN results. 


3. Make the integration step size used by Eq.(2.3) a function of the wave number 
to determine if improved accuracy can be obtained. 
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APPENDIX. TABULATED VALUES FOR 17(y), k2(y), and C(y) 
FOR DIFFERENT LAUNCH ANGLES AND RECEIVER DEPTHS. 


FREQUENCY = 50 HZ 
SOURCE DEPTH = 150 M 
n’(y) WITH POSITIVE GRADIENT 


RECEIVER DEPTH 


-0.0005393 
1.6446 


0.98673 
-0.0003760 


0.0016844 
-5.1370 





Table 18. Values of n?(y), k2(y), C(y) for 
different launch angles By and receiver depths y. 








FREQUENCY = 250 HZ 
SOURCE DEPTH = 150 M 
n’(y) WITH POSITIVE GRADIENT 


RECEIVER DEPTH 


0.0054445 


0.0054445 


0.0048032 





Table 19. Values of n?(y), k2(y), C(y) for 
different launch angles By and receiver depths y. 
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FREQUENCY = 1000 HZ 
SOURCE DEPTH = 150 M 
n’(y) WITH POSITIVE GRADIENT 


RECEIVER DEPTH 





Table 20. Values of n?(y), k2(y), C(y) for 
different launch angles Bo and receiver depths y. 















FREQUENCY = 50 HZ 
SOURCE DEPTH = 150 M 
n’(y) WITH NEGATIVE GRADIENT 


RECEIVER DEPTH 
SS 


1.0067 
0.00075794 
-2.3562 















1.0067 
0.00049734 
-1.5474 






0.00078545 
2.4438 














1.0134 
0.00062858 


1.0067 
0.00034048 





-0.0005238 








0.00005237 
-0.16295 


9.00091669 










1.00000 
0.00020924 
-0.65101 







0.00046983 
-1.4618 











1.00000 
0.00083289 
-2.5914 







Table 21. Values of n(v), k2(y), C(y) for 
different launch angles Bo and receiver depths y. 
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FREQUENCY = 250 HZ 
SOURCE DEPTH = 150 M 
n’(y) WITH NEGATIVE GRADIENT 


RECEIVER DEPTH 





Table 22. Values of n?(y), k2(y), G(y) for 
different launch angles Bo and receiver depths y. 















FREQUENCY = 1000 HZ 
SOURCE DEPTH = 150 M 
n’(y) WITH NEGATIVE GRADIENT 


RECEIVER DEPTH 


y=0 y=50 y= 100 y= 150 y=200 y= 250 y = 300 
meters Ineters meters meters meters meters meters 





Table 23. Values of n(y), k2(y), C(y) for 
different launch angles By and receiver depths y. 


58 





BETAO = 84 BETAO = 86 BETAO = 68 
i) 


=x 100 =X 100 = 100 
tt 200 4 200 200 


300 300 300 
6 $6 6 


o = 66 0 66 


0 66 


ZETA ZETA ZETA 
BETAO = 92 BETAO = 94 gpelnuee 


0 
= 100 =x 100 =X 100 
& im 
& 200 8 200 & 200 
300 900 300 
§ 0 § § 0 § 6 
ZETA ZETA ZETA 
FREQUENCY = §0 hertz 
SOURCE DEPTH = 160 meters 
OCEAN DEPTH = 300 meters 


SOUND SPEED, OCEAN SURFACE = 100 m's 
SOUND SPEED, SOURCE DEPTH = 1486 ms 





Figure 12. Plot of C(y) versus depth y for different launch angles Bo: 
n*(y) with positive gradient. 


P scalhied = 90 


=x 100 
200 


300 
-10 0 10 


ZETA 
BETAS = 98 
0 
=x 100 
& 
4 200 


300 
20 0 20 -20 6 20 
ZETA ZETA 


FREQUENCY = 260 hertz 

SOURCE DEPTH = 150 meters 

OCEAN DEPTH = 300 meters 

SOUND SPEED, OCEAN SURFACE = 1500 mis 
SOUND SPEED, SOURCE DEPTH = 1485 més 





Figure 13. Plot of C(y) versus depth y for different launch angles Bo: 
n’(y) with positive gradient. 








FREQUENCY = 1000 hertz 

SOURCE DEPTH = 150 meters 

OCEAN DEPTH = 900 meters 

SOUND SPEED, OCEAN SURFACE = 1500 ms 
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Figure 14. Plot of C(y) versus depth y for different launch angles Bo: 
n’(y) with positive gradient. 
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Figure 15. Plot of C(v) versus depth y for different launch angles Bo: 
n’(y) with negative gradient. 
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Figure 16. Plot of C(y) versus depth y for different launch angles Bo: 
n*(y) with negative gradient. 
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Figure 17. Plot of C(v) versus depth y for different launch angles Bo: 
n*(y) with negative gradient. 
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